library(lme4)
library(ggplot2)
library(ggeffects)
library(patchwork)
library(stargazer)
library(lavaan)
library(tidySEM)
library(dplyr)

load("wkdf_rep.rda")


# Trust #-----------------------------------------------------------------------------------------#
#
# Main analysis

f0 <- lmer(trust_sum2 ~ type2 +  (1 | cty.x),
           data = wk.df)
summary(f0)

f1 <- lmer(trust_sum2 ~ type2 +  (1 + type2 | cty.x),
           data = wk.df)
summary(f1)

f2 <- lmer(trust_sum2 ~ type2 +  sex + edu + (1 + type2 | cty.x),
           data = wk.df)
summary(f2)

f3 <- lmer(trust_sum2 ~ type2 +  sex + edu + age_g + as.factor(year_f) + (1 + type2 | cty.x),
           data = wk.df)
summary(f3)

f4 <- lmer(trust_sum2 ~ type2 + sex + edu + age_g + as.factor(year_f) + 
             log(gdppc) + democracy + postcomm + state_support + (1 + type2 | cty.x),
           data = wk.df)
summary(f4)

f5 <- lmer((trust1/5) ~ type2 + sex + edu + age_g + as.factor(year_f) + 
             log(gdppc) + democracy + postcomm + state_support + (1 + type2 | cty.x),
           data = wk.df)
summary(f5)

f7 <- lmer((trust3/5) ~ type2 + sex + edu + age_g + as.factor(year_f) + 
             log(gdppc) + democracy + postcomm + state_support + (1 + type2 | cty.x),
           data = wk.df)
summary(f7)


# Table 1
stargazer(f1, f2, f3, f4, f5, f7)


# Trust #-----------------------------------------------------------------------------------------#
#
# Heterodox without healing

f1a <- lmer(trust_sum2 ~ type2_3 + (1 + type2_3 | cty.x),
            data = wk.df)
summary(f1a)

f2a <- lmer(trust_sum2 ~ type2_3 + sex + edu + (1 + type2_3 | cty.x),
            data = wk.df)
summary(f2a)

f3a <- lmer(trust_sum5 ~ type2_3 + sex + edu + age_g + as.factor(year_f) + (1 + type2_3 | cty.x),
            data = wk.df)
summary(f3a)

f4a <- lmer(trust_sum2 ~ type2_3 + sex + edu + age_g + as.factor(year_f) + 
              log(gdppc) + democracy + postcomm + state_support  + (1 + type2_3 | cty.x),
            data = wk.df)
summary(f4a)

f5a <- lmer((trust1/5) ~ type2_3 + sex + edu + age_g + as.factor(year_f) + 
              log(gdppc) + democracy + postcomm + state_support + (1 + type2_3 | cty.x),
            data = wk.df)
summary(f5a)

f7a <- lmer((trust3/5) ~ type2_3 + sex + edu + age_g + as.factor(year_f) + 
              log(gdppc) + democracy + postcomm + state_support + (1 + type2_3 | cty.x),
            data = wk.df)
summary(f7a)


# Table E1
stargazer(f1a, f2a, f3a, f4a, f5a, f7a)


# Trust #-----------------------------------------------------------------------------------------#
#
# Heterodox of four categories

wk.df$type <- factor(wk.df$type, levels = c("Relig_ex", "Relig_in", "Rational", "Heterodox"))
wk.df$type <- relevel(wk.df$type, ref = "Relig_in")

a1 <- lmer(trust_sum2 ~ type +  (1 + type | cty.x),
           data = wk.df)
summary(a1)

a4 <- lmer(trust_sum2 ~ type + sex + edu + age_g + as.factor(year_f) + 
             log(gdppc) + democracy + postcomm + state_support + (1 + type | cty.x),
           data = wk.df)
summary(a4)

a5 <- lmer((trust1/5) ~ type + sex + edu + age_g + as.factor(year_f) + 
             log(gdppc) + democracy + postcomm + state_support + (1 + type | cty.x),
           data = wk.df)
summary(a5)

a7 <- lmer((trust3/5) ~ type + sex + edu + age_g + as.factor(year_f) + 
             log(gdppc) + democracy + postcomm + state_support + (1 + type | cty.x),
           data = wk.df)
summary(a7)


# Table F1
stargazer(a1, a4, a5, a7)


# Anti-Science #----------------------------------------------------------------------------------#

wk.df$type2 <- relevel(as.factor(wk.df$type2), ref = "Rational")

f1 <- lmer(sc ~ type2 + (1 + type2 | cty.x),
           data = wk.df)
f1a <- lmer(sc1 ~ type2 + (1 + type2 | cty.x),
            data = wk.df)
f1b <- lmer(sc2 ~ type2 + (1 + type2 | cty.x),
            data = wk.df)
summary(f1); summary(f1a); summary(f1b)

f2 <- lmer(sc ~ type2 + sex + edu + (1 + type2 | cty.x),
           data = wk.df)
summary(f2)

f3 <- lmer(sc ~ type2 +  sex + edu + age_g + year_f + (1 + type2 | cty.x),
           data = wk.df)
summary(f3)

f4 <- lmer(sc ~ type2 + sex + edu + age_g + as.factor(year_f) + 
             log(gdppc) + democracy + postcomm + state_support + (1 + type2 | cty.x),
           data = wk.df)
f4a <- lmer(sc1 ~ type2 + sex + edu + age_g + as.factor(year_f) + 
              log(gdppc) + democracy + postcomm + state_support + (1 + type2 | cty.x),
            data = wk.df)
f4b <- lmer(sc2 ~ type2 + sex + edu + age_g + as.factor(year_f) + 
              log(gdppc) + democracy + postcomm + state_support + (1 + type2 | cty_iso3),
            data = wk.df)
summary(f4); summary(f4a); summary(f4b)


# Mediation Analysis #----------------------------------------------------------------------------#
wk.df$heterodox <- ifelse(wk.df$type2 == "Heterodox", 1, 0)
wk.df$rational <- ifelse(wk.df$type2 == "Rational", 1, 0)
wk.df$religious <- ifelse(wk.df$type2 == "Religious", 1, 0)

m1 <- '
  # mediator
  sc ~ a1*religious + a2*rational

  # outcome
  trust1 ~ b*sc + b1*religious + b2*rational

  # mediated effect
  me_rel := a1*b
  me_rat := a2*b

  # direct effect
  de_rel := b1
  te_rat := b2
  
  # total effect
  te_rel := b1 + a1*b
  te_rat := b2 + a2*b
'

f1 <- sem(m1, data = wk.df)
summary(f1)

l <- matrix(NA, 3, 4); l
l[1,1] <- "religious"
l[3,1] <- "rational"
l[2,4] <- "trust1"
l[2,2] <- "sc"
l

newlab <- c("rational", "religious", "anti-science", "trust in leg.")
p <- prepare_graph(f1,  layout = l, 
                   angle = 90, variance_diameter = .5)
nodes(p)$label <- newlab  

p %>%
  color_pos_edges("green") %>%
  color_neg_edges("red") %>%
  color_var("black") %>%
  alpha_var(.3) %>%
  plot()


m2 <- '
  # mediator
  sc ~ a1*religious + a2*rational

  # outcome
  trust3 ~ b*sc + b1*religious + b2*rational

  # mediated effect
  me_rel := a1*b
  me_rat := a2*b

  # direct effect
  de_rel := b1
  de_rat := b2
  
  # total effect
  te_rel := b1 + a1*b
  te_rat := b2 + a2*b
'

f2 <- sem(m2, data = wk.df)
summary(f2)

l <- matrix(NA, 3, 4); l
l[1,1] <- "religious"
l[3,1] <- "rational"
l[2,4] <- "trust3"
l[2,2] <- "sc"
l

newlab <- c("rational", "religious", "anti-science", "trust in chu.")
p <- prepare_graph(f2,  layout = l, 
                   angle = 90, variance_diameter = .5)
nodes(p)$label <- newlab  

p %>%
  color_pos_edges("green") %>%
  color_neg_edges("red") %>%
  color_var("black") %>%
  alpha_var(.3) %>%
  plot()
# 4*5.5

